####################################
# Code to Replicate Charts Figures and Tables
# Trust in Public Policy Algorithms
# Journal of Politics
# Criminal Justice Experiments
# Study 3
# Ryan Kennedy, Philip D. Waggoner, and Matthew Ward
# 2021
####################################

#### Load needed libraries
library(tidyverse)
library(here)
library(arm)
library(stargazer)

#### Load the data
dta <- read_csv(here("HybridExperiment.csv"))

#### Code the treatment
rdta <- dta %>%
  rowwise() %>%
  mutate(s1outcome = sum(`s1 agree high`, `s1 agree low`, 
                         `s1 disagree high`, `s1 disagree low`, na.rm = T),
         s2outcome = sum(`s2 agree high`, `s2 agree low`, 
                         `s2 disagree high`, `s2 disagree low`, na.rm = T),
         s3outcome = sum(`s3 agree high`, `s3 agree low`, 
                         `s3 disagree high`, `s3 disagree low`, na.rm = T),
         s4outcome = sum(`s4 agree high`, `s4 agree low`, 
                         `s4 disagree high`, `s4 disagree low`, na.rm = T),
         s5outcome = sum(`s5 agree high`, `s5 agree low`, 
                         `s5 disagree high`, `s5 disagree low`, na.rm = T),
         s6outcome = sum(`s6 agree high`, `s6 agree low`, 
                         `s6 disagree high`, `s6 disagree low`, na.rm = T),
         s7outcome = sum(`s7 agree high`, `s7 agree low`, 
                         `s7 disagree high`, `s7 disagree low`, na.rm = T)) %>%
  ungroup() %>%
  gather(scenario, outcome, s1outcome:s7outcome) %>%
  mutate(outcome = ifelse(outcome == 0, NA, outcome)) %>%
  mutate(treatment = case_when(scenario == "s1outcome" & !is.na(`s1 agree high`) ~ "Agree High",
                               scenario == "s1outcome" & !is.na(`s1 agree low`) ~ "Agree Low",
                               scenario == "s1outcome" & !is.na(`s1 disagree high`) ~ "Disagree High",
                               scenario == "s1outcome" & !is.na(`s1 disagree low`) ~ "Disagree Low",
                               scenario == "s2outcome" & !is.na(`s2 agree high`) ~ "Agree High",
                               scenario == "s2outcome" & !is.na(`s2 agree low`) ~ "Agree Low",
                               scenario == "s2outcome" & !is.na(`s2 disagree high`) ~ "Disagree High",
                               scenario == "s2outcome" & !is.na(`s2 disagree low`) ~ "Disagree Low",
                               scenario == "s3outcome" & !is.na(`s3 agree high`) ~ "Agree High",
                               scenario == "s3outcome" & !is.na(`s3 agree low`) ~ "Agree Low",
                               scenario == "s3outcome" & !is.na(`s3 disagree high`) ~ "Disagree High",
                               scenario == "s3outcome" & !is.na(`s3 disagree low`) ~ "Disagree Low",
                               scenario == "s4outcome" & !is.na(`s4 agree high`) ~ "Agree High",
                               scenario == "s4outcome" & !is.na(`s4 agree low`) ~ "Agree Low",
                               scenario == "s4outcome" & !is.na(`s4 disagree high`) ~ "Disagree High",
                               scenario == "s4outcome" & !is.na(`s4 disagree low`) ~ "Disagree Low",
                               scenario == "s5outcome" & !is.na(`s5 agree high`) ~ "Agree High",
                               scenario == "s5outcome" & !is.na(`s5 agree low`) ~ "Agree Low",
                               scenario == "s5outcome" & !is.na(`s5 disagree high`) ~ "Disagree High",
                               scenario == "s5outcome" & !is.na(`s5 disagree low`) ~ "Disagree Low",
                               scenario == "s6outcome" & !is.na(`s6 agree high`) ~ "Agree High",
                               scenario == "s6outcome" & !is.na(`s6 agree low`) ~ "Agree Low",
                               scenario == "s6outcome" & !is.na(`s6 disagree high`) ~ "Disagree High",
                               scenario == "s6outcome" & !is.na(`s6 disagree low`) ~ "Disagree Low",
                               scenario == "s7outcome" & !is.na(`s7 agree high`) ~ "Agree High",
                               scenario == "s7outcome" & !is.na(`s7 agree low`) ~ "Agree Low",
                               scenario == "s7outcome" & !is.na(`s7 disagree high`) ~ "Disagree High",
                               scenario == "s7outcome" & !is.na(`s7 disagree low`) ~ "Disagree Low")) %>%
  mutate(treatment = factor(treatment, levels = c("Agree Low", "Agree High", "Disagree High", "Disagree Low")))


#### Code the other IVs
rdta <- rdta %>%
  mutate(tia = (trustinautomation_1 + trustinautomation_2 + trustinautomation_3 +
                  (8 - trustinautomation_4) + (8 - trustinautomation_5) + 
                  trustinautomation_6 + trustinautomation_7 + 
                  (8 - trustinautomation_8))/8,
         female = ifelse(gender == 2, 1, 0),
         edvalue = case_when(education == 1 ~ 1,
                             education == 2 ~ 2,
                             education == 3 ~ 2,
                             education == 4 ~ 3,
                             education == 5 ~ 3,
                             education == 6 ~ 3,
                             education == 7 ~ 4,
                             education == 8 ~ 4),
         partisanship = case_when(political_party == 7 | political_party == 4 ~ 4,
                                  political_party == 6 | political_party == 3 ~ 3,
                                  political_party == 8 | political_party == 5 ~ 5,
                                  political_party == 9 ~ 6,
                                  political_party == 2 ~ 2,
                                  political_party == 10 ~ 7,
                                  political_party == 1 ~ 1))

#### Table of summary statistics
# Table A4 in SI
srdta <- rdta %>%
  filter(Finished == 1) %>%
  dplyr::select(outcome, treatment, tia, age, edvalue, female, partisanship) %>%
  mutate(agreeLow = ifelse(treatment == "Agree Low", 1, 0),
         agreeHigh = ifelse(treatment == "Agree High", 1, 0),
         dHigh = ifelse(treatment == "Disagree High", 1, 0),
         dLow = ifelse(treatment == "Disagree Low", 1, 0)) %>%
  dplyr::select(outcome, agreeHigh, agreeLow, dHigh, dLow, tia, age, edvalue, female, partisanship)

srdta <- data.frame(srdta)

sink(here("hybridSummaries.tex"))
stargazer(srdta, style="ajps", 
          title="Summary Statistics", 
          covariate.labels=c("Outcome", "Algorithm & Judge High", 
                             "Algorithm & Judge Low", "Algorithm High, Judge Low",
                             "Algorithm Low, Judge High",
                             "Trust in Automation", "Age",
                             "Education", "Female",
                             "Partisanship")
)
sink()


##### Baseline models of treatment conditions on outcomes
# Figure 3 in main text and Table A10 in SI
model1 <- lmer(outcome ~ treatment + (1|ResponseId) + (1|scenario), 
               data = rdta, control = lmerControl(optimizer = "Nelder_Mead"))
summary(model1)
# Table A18 in SI
model2 <- lmer(outcome ~ treatment + tia + age + edvalue + female + partisanship 
               + (1|ResponseId) + (1|scenario), 
               data = rdta, control = lmerControl(optimizer = "Nelder_Mead"))
summary(model2)


# Create plot of model 1 results
# Figure 3 in main text -- left side
mod1sum <- summary(model1)
mod1sim <- sim(model1, n.sims = 1000)
mod1simcoef <- as_tibble(data.frame(coef(mod1sim)$fixef))
mod1simcoef$simulation <- rownames(mod1simcoef)
mod1simcoeflong <- gather(mod1simcoef, condition, coefficient, X.Intercept.:treatmentDisagree.Low)
mod1simcoeflong$condition[mod1simcoeflong$condition == "treatmentAgree.High"] <- "Algorithm & Judge High"
mod1simcoeflong$condition[mod1simcoeflong$condition == "treatmentDisagree.High"] <- "Algorithm High, Judge Low"
mod1simcoeflong$condition[mod1simcoeflong$condition == "treatmentDisagree.Low"] <- "Algorithm Low, Judge High"
mod1simcoeflong$condition <- factor(mod1simcoeflong$condition, levels = c("X.Intercept.", "Algorithm & Judge High", "Algorithm Low, Judge High", "Algorithm High, Judge Low"))

m1b95 <- mod1simcoeflong %>%
  filter(condition == "Algorithm & Judge High" | condition == "Algorithm High, Judge Low" | condition == "Algorithm Low, Judge High") %>%
  group_by(condition) %>%
  summarize(mcond = mean(coefficient, na.rm = T),
            sdcond = sd(coefficient, na.rm = T)) %>%
  mutate(hi = mcond + 1.96*sdcond,
         lo = mcond - 1.96*sdcond) %>% 
  ggplot() +
  geom_point(aes(x = condition, y = mcond), size = 3, color = "darkgreen") +
  geom_errorbar(aes(x = condition, ymin = lo, ymax = hi), color = "darkgreen") +
  geom_hline(aes(yintercept = 0)) +
  labs(x = "Condition", y = "ATE", title = "Likelihood of Re-offense") +
  coord_flip() + theme_bw()
m1b95

# Simulation of linear model
# Figure 3 in main text -- right side
mlin95 <- mod1simcoeflong %>%
  filter(condition == "Algorithm & Judge High" | condition == "Algorithm High, Judge Low" | condition == "Algorithm Low, Judge High") %>%
  group_by(condition) %>%
  summarize(mcond = mean(coefficient, na.rm = T),
            sdcond = sd(coefficient, na.rm = T)) %>%
  mutate(hi = mcond + 1.96*sdcond,
         lo = mcond - 1.96*sdcond) %>% 
  ungroup() %>%
  mutate(mlin = case_when(condition == "Algorithm & Judge High" ~ 0.9,
                          condition == "Algorithm High, Judge Low" ~ 0.3,
                          condition == "Algorithm Low, Judge High" ~ 0.6),
         hilin = mlin + (1.96 * 0.065),
         lolin = mlin -(1.96 * 0.065)) %>%
  ggplot() +
  geom_point(aes(x = condition, y = mcond, color = "Results"), size = 3) +
  geom_errorbar(aes(x = condition, ymin = lo, ymax = hi, color = "Results")) +
  geom_point(aes(x = condition, y = mlin, color = "Linear Model"), size = 3) +
  geom_errorbar(aes(x = condition, ymin = lolin, ymax = hilin, color = "Linear Model")) +
  geom_hline(aes(yintercept = 0)) +
  scale_color_manual(values = c("dark green", "red")) +
  labs(x = "Condition", y = "ATE", title = "Likelihood of Re-offense", color = "Model") +
  coord_flip() + theme_bw() + theme(legend.position = c(0.8, 0.85))
mlin95

# Formatting and saving Table A10 in SI
insertrow <- function(existingDF, newrow, r) {
  existingDF[seq(r + 1,nrow(existingDF) + 1),] <- existingDF[seq(r,nrow(existingDF)),]
  existingDF[r,] <- newrow
  existingDF
}
# Create a baseline table
Tables <- stargazer(model1, style = "ajps", 
                    title = "Models for Likelihood of Committing a Crime if Released", 
                    dep.var.labels.include = FALSE, 
                    covariate.labels = c("Algorithm & Judge High", "Algorithm High, Judge Low", 
                                         "Algorithm Low, Judge High")
)
# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)
# Find where you want to put in the random effect. In our case, this is right after the last fixed effect.
r <- 16
# Create some standard label lines
randomeffect <- "{\\bf Random Effect} & \\\\"
hline <- "\\hline"
newline <- "\\\\"
# Now insert those label lines where they are needed
Tables <- insertrow(Tables, hline, r)
Tables <- insertrow(Tables,randomeffect,r + 1)
Tables <- insertrow(Tables,hline,r + 2)
# Get number of unique values in each grouping
num.participants1 <- sapply(ranef(model1),nrow)[1]
num.scenarios1 <- sapply(ranef(model1),nrow)[2]
# Get standard deviation of the random effect
stddev.model1.participants <- attributes(VarCorr(model1)$ResponseId)$stddev
stddev.model1.scenarios <- attributes(VarCorr(model1)$scenario)$stddev
# Create a LaTex character row for the random effect
number.of.participants <- paste("\\# of Participants & ", num.participants1, "\\\\")
stddev.participants <- paste("Participant Standard Deviation & ", round(stddev.model1.participants, 3), "\\\\")
number.of.scenarios <- paste("\\# of Scenarios & ", num.scenarios1, "\\\\")
stddev.scenarios <- paste("Scenario Standard Deviation & ", round(stddev.model1.scenarios, 3), "\\\\")
# Add these lines to the table
Tables <- insertrow(Tables,number.of.participants,r + 3)
Tables <- insertrow(Tables,stddev.participants,r + 4)
Tables <- insertrow(Tables,newline,r + 5)
Tables <- insertrow(Tables,number.of.scenarios,r + 6)
Tables <- insertrow(Tables,stddev.scenarios,r + 7)
Tables <- insertrow(Tables, hline, r + 8)
# Write the table to a file that can be inserted into the document
write.table(Tables,file = here("REtable(Exp5).tex"),
            sep = "",row.names = FALSE,na = "", quote = FALSE, col.names = FALSE)

# Formatting and saving Table A18 in SI
Tables <- stargazer(model2, style = "ajps", 
                    title = "Models for Likelihood of Committing a Crime if Released", 
                    dep.var.labels.include = FALSE, 
                    covariate.labels = c("Algorithm & Judge High", "Algorithm High, Judge Low", 
                                         "Algorithm Low, Judge High", "Trust in Automation", 
                                         "Age","Education", "Female", "Republican Partisanship")
)
# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)
# Find where you want to put in the random effect. In our case, this is right after the last fixed effect.
r <- 27
# Create some standard label lines
randomeffect <- "{\\bf Random Effect} & \\\\"
hline <- "\\hline"
newline <- "\\\\"
# Now insert those label lines where they are needed
Tables <- insertrow(Tables, hline, r)
Tables <- insertrow(Tables,randomeffect,r + 1)
Tables <- insertrow(Tables,hline,r + 2)
# Get number of unique values in each grouping
num.participants2 <- sapply(ranef(model2),nrow)[1]
num.scenarios2 <- sapply(ranef(model2),nrow)[2]
# Get standard deviation of the random effect
stddev.model2.participants <- attributes(VarCorr(model2)$ResponseId)$stddev
stddev.model2.scenarios <- attributes(VarCorr(model2)$scenario)$stddev
# Create a LaTex character row for the random effect
number.of.participants <- paste("\\# of Participants & ", num.participants2, "\\\\")
stddev.participants <- paste("Participant Standard Deviation & ", round(stddev.model2.participants, 3), "\\\\")
number.of.scenarios <- paste("\\# of Scenarios & ", num.scenarios2, "\\\\")
stddev.scenarios <- paste("Scenario Standard Deviation & ", round(stddev.model2.scenarios, 3), "\\\\")
# Add these lines to the table
Tables <- insertrow(Tables,number.of.participants,r + 3)
Tables <- insertrow(Tables,stddev.participants,r + 4)
Tables <- insertrow(Tables,newline,r + 5)
Tables <- insertrow(Tables,number.of.scenarios,r + 6)
Tables <- insertrow(Tables,stddev.scenarios,r + 7)
Tables <- insertrow(Tables, hline, r + 8)
# Write the table to a file that can be inserted into the document
write.table(Tables,file = here("REtable(Exp5)wcov.tex"),
            sep = "",row.names = FALSE,na = "", quote = FALSE, col.names = FALSE)


#### Create balance table
# Table A7 in SI
rdta <- rdta %>%
  mutate(alghijdglo = ifelse(treatment == "Disagree High", 1, 0),
         alglojdghi = ifelse(treatment == "Disagree Low", 1, 0),
         alghijdghi = ifelse(treatment == "Agree High", 1, 0))

bal1 <- lm(alghijdglo ~ age + edvalue + female + partisanship, 
           data = rdta)
summary(bal1)

bal2 <- lm(alglojdghi ~ age + edvalue + female + partisanship, 
           data = rdta)
summary(bal2)

bal3 <- lm(alghijdghi ~ age + edvalue + female + partisanship, 
           data = rdta)
summary(bal3)

# Create and save Table A7
# Create a baseline table
Tables <- stargazer(bal1, bal2, bal3, style="ajps", 
                    title="Hybrid Algorithm Balance Table", 
                    dep.var.labels.include = FALSE, 
                    covariate.labels=c("Age","Education","Female", "Republican Partisanship")
)
# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)
# Write the table to a file that can be inserted into the document
write.table(Tables,file=here("ajBalTable.tex"),
            sep="",row.names= FALSE,na="", quote = FALSE, col.names = FALSE)


